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We present an algorithm for the simulation of the exact real-time dynamics of classical many-body 
systems with discrete energy levels. In the same spirit of kinetic Monte Carlo methods, a stochastic 
solution of the master equation is found, with no need to dehne any other phase-space construction. 
However, unlike existing methods, the present algorithm does not assume any particular statistical 
distribution to perform moves or to advance the time, and thus is a unique tool for the numerical 
exploration of fast and ultra-fast dynamical regimes. By decomposing the problem in a set of two- 
level subsystems, we find a natural variable step size, that is well defined from the normalization 
condition of the transition probabilities between the levels. We successfully test the algorithm with 
known exact solutions for non-equilibrium dynamics and equilibrium thermodynamical properties of 
Ising-spin models in one and two dimensions, and compare to standard implementations of kinetic 
Monte Carlo methods. The present algorithm is directly applicable to the study of the real time 
dynamics of a large class of classical markovian chains, and particularly to short-time situations 
where the exact evolution is relevant. 

PACS numbers: 05.10.Ln, 02.50.Ga, 05.70.Ln 


INTRODUCTION 

Many body classical models with discrete energy lev¬ 
els, such as Ising-spin systems, are particular examples 
of markovian chains [ij, whose growing interest includes 
fields as diverse as condensed matter physics, biology 
0 and economics Q. Despite intense research, exact 
results for these systems are rare in statistical physics, 
even for the most simple Hamiltonians Q. In this con¬ 
text, Monte Carlo (MC) numerical calculations are often 
considered as a fundamental benchmark for theories and 
experiments Q. 

While MC simulations usually provide accurate results 
for static properties of interacting discrete-variable mod¬ 
els, the situation is different regarding their dynamical 
evolution, which, lacking a first-principles equation of 
motion as in continuous-variable systems, should b e g en- 
erally described by a stochastic master equation [l|, [gj. 
The latter expresses the probability distribution 'P{X,t) 
of a given state X at time t, in the formQ 

^ W{X\Y)P{Y, t)-J2 W(Y\X)V(Y, t) (1) 

Y Y 

where W{Y\X) is the transition rate from state X to 
state Y, in units of inverse time. 

In model with discrete variables, where the states form 
a numerable set, the common requirement for a dynam¬ 
ical MC algorithm is to reach asymptotically the equi¬ 
librium state, where the master equation fulfills detailed 
balance Q . As a result Monte Carlo algorithms are usu¬ 
ally based on the equilibrium (e.g., time-independent) 
transition probabilities between states, instead of the 


time-varying probabilities resulting from the general so¬ 
lution of equation O- The standard Monte Carlo step 
(MCS) that is used as the time step in most algorithms 
thus measures just the extent of random exploration over 
the configuration space and has no direct relation with 
physical time. In general, this can result in significant 
deviations between the MC dynamics and the dynami¬ 
cal behavior described by the master equation. However 
some equilibrium algorithms are known to reproduce suc¬ 
cessfully certain dynamical laws. For example, this is 
the case of the Metropolis algorithm that predicts the 
m ~ scaling for the magnetization m of the 2D Ising 
model after a subcritical quench Q. 

So far, the most important bridge between MCS and 
physical time has been built by a class of algorithms usu¬ 
ally called dynamic or kinetic Monte Carlo (KMC) [^, . 

KMC algorithms use the information about the transi¬ 
tion rates W(X\Y) to select the new updates, thus as¬ 
signing to this process a real time related to the inverse 
rates. 

More specifically, KMC algorithms use the fact that 
the average time between two consecutive events in the 
system is of the order of (At) ^ R(X)~^, where R{X) = 
W(Y\X) is the total sum of all rates of individual 
processes the system can undergo from a given state 
X 0. Therefore, single time step is updated in a real¬ 
istic way using a Poissonian distribution, by the expres¬ 
sion At = —i?(Ar)“^log(a;), with x being a uniformly dis¬ 
tributed variable between 0 and I. This trick allows one 
to map the simulation steps with a real time that is phys¬ 
ically meaningful, and has become the current standard 
for numerical calculations of the dynamics of discrete- 
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variable models. The KMC step is then completed by the 
execution of the process that has been selected following 
a specific rule. The choice of this specihc rule have pro¬ 
duced different KMC schemes: the so called first-reaction 
method for example, selects always the process with 
the fastest rate, while in the most commonly used BKL 
or Gillespie algorithm [13 the probability of selecting a 
process is a linear function of the rates. 

As can be inferred from the discussion above, all stan¬ 
dard KMC methods follow a Markov chain kinetics, sam¬ 
pling correctly from the (usually unknown) solution of 
the master equation, and so producing stochastic trajec¬ 
tories along the actual time axis. These single trajec¬ 
tories, however, are very accurate as far as time scales 
remain larger than (At). At times of the order of con¬ 
secutive events, trajectories are not expected to repro¬ 
duce the exact solution of the master equation in the 
time axis. This loss of accuracy at small times prevents, 
for example, the inclusion in the KMC dynamics of any 
time-dependent parameter whose variation is of the order 
of (At). A reliable numerical technique capable of repro¬ 
ducing the master equation kinetics for fast and ultra-fast 
regimes is still lacking. 

In this work we present a new algorithm for addressing 
the latter problem, that is based on the numerical solu¬ 
tion of the master equation. The main requirements are 
that (i) the system can be decomposed into a set of N 
two-level subsystems, and (ii) any dynamical evolution is 
realized by sequential transitions within these individual 
subsystems. In the following we refer to these N transi¬ 
tions as minimal processes. Condition (i) is the standard 
form of any Hamiltonians with Ising-like spins, however, 
it can be also made to apply to, e.g., classical mixtures 
on lattices or any generic Potts models. Condition (ii) is 
equivalent to the well known single-spin-flip update pro¬ 
cedure, which is widely used for dynamical calculations of 
discrete models. As we show below, conditions (i) and (ii) 
can be fulfilled in any model with discrete-variables. We 
expect that the algorithm will be of particular value for, 
e.g., short-time critical dynamics of interacting classical 
models 111, phase order kinetics [l^ and driven systems 
in oscillating fields [13. However, its validity is not re¬ 
stricted to physical systems, nor to short times, and may 
in principle be used as an alternative to Metropolis or 
KMC simulations in a large number of markovian chains 
of different nature. 


EVENT-DRIVEN ALGORITHM 

Without loss of generality, in the following we present 
the algorithm in terms of Ising spins. The idea behind the 
scheme is the following: given an initial configuration of 
the interacting spins, within the characteristic time scale 
Ted (which is a priori unknown) associated with the flip 
of a single spin from that specific configuration, the time 


evolution of the whole system is described by a set of N 
independent reduced master equations. By solving the 
latter, the exact time dependent probability Pi{At) for 
spin i to flip is obtained analytically for each spin i at 
any time At, where At is the time interval since the pre¬ 
vious spin flip. In turn, the condition = 1 

defines the value of ted consistent with the single spin 
flip for that given configuration. Once ted is defined, the 
algorithm proceeds with evaluating all Pi (At) at time 
At = TED and uses them to update the configuration. 
This concludes a step of the algorithm. The whole pro¬ 
cedure is then repeated. 

At the beginning of each step, we consider each spin 
i occupying level o while level / is initially free, so that 
= 0) = 1 and P/ (0) = 0, where P°^^{t) corre¬ 
sponds to the occupation probabilities of the two levels. 
Within Ted , these occupation probabilities fulfill the rate 
equations 


= r{‘>p/{t)-rfpf{t) 

= Tfpnt) - T^P/ it) ( 2 ) 

where and are the transition rates of the two-level 
subsystem of spin i, depending on the energy value 
of the levels and the physical nature of the system. The 
limit of infinite time corresponds to Boltzmann occupa¬ 
tion probabilities P°’^(oo) = fz, where Z is the 

partition function of the two-level subsystem, j3~^ = ksT 
and fcs is the Boltzmann constant. 

With these conditions the transition probability 
Pi{At), i.e., from o to /, can be written in the form 


dprit) 

dt 

dPRt) 

dt 


P,{At) = p/{At) 


P/{oo) 


1 _ g-(r°^-erf°)At 


( 3 ) 


where P/(oo) = ^]. As usually done in lit¬ 

erature, in the following we assume that the characteris¬ 
tic frequency T = r°'^-|-r{° is constant in the system, and 
1 /r is adopted as the unit of time [3 • From Eq. ([3]) and 
by applying the normalization condition Pi(TED) = 1 
given above, we obtain ted as 

Tted = -In [1 - P,-^] , (4) 

with P, = -P/(oo)- 

Each step of the algorithm starts with the calculation, 
for each spin i, of the energy difference AEi = e( — E° 
associated to flipping the spin. From this, the value of 
P* can be calculated. If P* > 1 then the value of ted for 
the current step becomes that of expression Eq. ([4]). If 
P* < 1, we choose to set ted = 1 (see below). Once ted 
is evaluated, the sites are updated with the correspond¬ 
ing probability Eq. ([3]) with At = ted, resulting in an 
average of one spin flip. Consequently, the total time of 
the simulation is now incremented by ted- 
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A value of P* less than or equal to one, means that 
the system will never reach a time for which, in average, 
one spin is flipped. This is the well-known situation in 
which finite systems freeze, and the dynamics arrests, af¬ 
ter reaching a stable configuration at sufficiently low tem¬ 
peratures [one example is shown below when discussing 
Fig. [2] (inset)]. In most cases, this condition should sug¬ 
gest the end of the calculation, since the system will never 
evolve after reaching this state. However, for problems in 
which the energy can change independently of the config¬ 
uration (e.g., time-dependent Hamiltonians), this freez¬ 
ing could be temporary and, consequently, ted should be 
set to a constant value when P* < 1. The value ted = 1 
is just a conventional number, since it has to be tailored 
to well-capture the time scale associated to the energy 
changes in the problem at hand. 

We note that ted corresponds to a discretization of real 
time and in general varies from step to step, as it is linked 
to the elementary changes of the system. In turn, the 
latter depends only on the microscopic interactions in the 
Hamiltonian and the specific spin configuration. Thus, 
since the whole algorithm directly deals with the exact 
real time, when conditions (i) and (ii) above are satisfied 
we expect that the results of the numerical simulation 
will reproduce well those of the exact master equation at 
all time scales. 

The role of ted can also be seen as a coarsening of the 
dynamics to the next physically meaningful time value, 
calculated exactly, and not generated from a distribution 
function as in KMC schemes. This time coarsening rep¬ 
resents the stochastic counterpart of that in event-driven 
molecular dynamics approaches [1^, thus corresponding 
to the waiting time connecting two consecutive (stochas¬ 
tic) events. Therefore we refer to the above-described 
scheme as the Event-Driven algorithm (ED). This algo¬ 
rithm is composed of two serial loops of size N, firstly 
performing the calculation of P* and secondly updating 
the minimal processes with the corresponding probabil¬ 
ity. Consequently the ED step is of complexity 0{N), 
that is, it scales linearly with the number of two-level 
subsystems, which is the same as the Metropolis Monte 
Carlo step. 

In general, for any discrete-variable markovian chain, 
starting from a given configuration, the dynamical step 
is a rule selecting the next configuration among N possi¬ 
ble choices. In terms of the ED scheme, the latter means 
that for any markovian chain one can build the set of N 
two-levels equations. The only input of the algorithm is 
the list of the N energy differences corresponding to each 
one of the possible choices of configurations. The main 
idea of the ED scheme relies on the very commonly used 
approximation that many coupled equations can be de¬ 
coupled for the very short time scale in which the system 
performs what we call a minimal process. Using this fact 
the algorithm finds the characteristic time ted for which 
only one minimal process is likely to happen. 


Consider, for instance, the g-levels Potts model, in 
which each spin can be in one of the q > 2 available 
states. For a system of Ng spins, this model will im¬ 
ply a number of subsystems of N = (q — l)Ns, since, for 
any given configuration, a minimal process consists in the 
transition of one of the Ng spins to one of its ((? — 1) avail¬ 
able states. Each of these N possible transitions is identi- 
Hed with a minimal process by the ED algorithm, though 
for this model those minimal processes corresponding to 
the same spin are excluded. Thus, we just need to eval¬ 
uate the energy difference associated to each of these N 
transitions. 

In general, the number of minimal processes is not even 
forced to be constant along consecutive steps, as is the 
case for example in the lattice gas model. The latter 
consists of particles that occupy certain positions in a 
lattice, and are able to move only to first-neighbouring 
empty sites. The minimal processes here should be taken 
as the set of all single possible moves that particles can 
perform. For a very diluted configuration, this number 
of subsystems is then N = ZNp, where Np and Z are 
the number of particles and the coordination of the lat¬ 
tice, respectively. However, when two particles become 
nearest neighbours, the number of minimal processes N 
is reduced. 


NUMERICAL TESTS 

In the following we implement the ED algorithm. Its 
accuracy is tested in a dynamical problem whose exact 
solution is known, and further compared to that of a 
state-of-the-art A-fold KMC algorithm, implemented via 
the KMCLib library 0. Further tests are also pre¬ 
sented to show the consistency of the ED scheme with 
well known equilibrium and dynamical behaviors of the 
2D Ising model, while discussing some specific features 
of the method. 


Glauber exact solution 

We start by testing the algorithm in the exploration 
of the temporal evolution of the local magnetization in a 
linear spin chain (ID) following a quench. Before general 
tests involving averages for the total magnetization, we 
compare here the predictions of the algorithm to the ex¬ 
act solution of the full master equation, as obtained by 
Glauber [^. Up to our knowledge, this remains the more 
complex discrete-variable statistical system for which the 
local magnetization dynamics has been analytically ob¬ 
tained, in the full range of time scales and for arbitrary 
initial conditions. In turn, this analysis for the non¬ 
equilibrium properties of the local order parameter is the 
most complex test to which the algorithm can be sub¬ 
jected. 
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FIG. 1: Temporal evolution of the average local magneti¬ 
zation {Si{t)) for ten different sites (from top to bottom 
i = 56, 57, 58,..., 65) of an Ising chain of L = 110 sites at 
T = 0.1. The initial condition was set to (Si(0)} = 1 for the 
ten central spins 51 < i < 60, and (S'i(O)) = 0 otherwise. 
Dots are the outcome of numerical simulations using Event- 
Driven algorithm (a), and KMC algorithm fitting the data for 
large (b) and small (c) spin index. Solid lines are the corre¬ 
sponding exact analytical solution obtained by R. J. Glauber 
in Ref. 



FIG. 2: Equilibrium magnetization m (blue) and energy 
E (red) obtained with the algorithm (triangles) and with 
Metropolis calculations (squares), by a slow annealing of the 
system at different temperatures. Units of energy and tem¬ 
perature are J and J/fcs, respectively, where J is the coupling 
constant of the Ising Hamiltonian. Magnetization is the av¬ 
erage value of the spins. The dashed line corresponds to the 
exact transition temperature of the infinite system. The inset 
is the temperature dependence of the average time step of the 
Event-Driven algorithm. 

all the sites at once. Moreover, as can be easily noticed 
from the figure, numerical and exact curves correspond¬ 
ing to sites near the edge of the block, are completely 
impossible to collapse by solely a rescaling of the time 
axis. 


Figure [TJi shows the time evolution of an Ising chain 
with L = 110 sites, where the initial state comprises 
a block of 10 parallel spins in the center, while the re¬ 
maining 100 are in a disordered state (see the caption 
for details). The figure shows a perfect agreement be¬ 
tween the exact analytical solution (continuous lines) and 
the numerical results from ED algorithm. Worth noting, 
this agreement occurs not only for the asymptotic, long¬ 
time regime, but also for very short times, where the 
system is strongly out of equilibrium and the functional 
dependence of the local magnetization on time is non¬ 
trivial. This confirms that the ED algorithm successfully 
accounts for the actual master-equation solution, accu¬ 
rately reproducing the trajectories in the real time axis, 
even for scales of the order of single flips. 

For comparison, figures [T)d and [TJ: shows the best fits 
for the outcome of the KMC algorithm in the same prob¬ 
lem. By adjusting the time scale with a free parameter, 
a reasonable fit can be found at short times for the local 
magnetization of sites far from (panel b), or deep into 
(panel c) the central ordered block of the initial chain 
configuration. While this rescaling is valid, it is impossi¬ 
ble to find a single rescaling parameter successfully fitting 


2D Isind model 

We now focus on equilibrium properties. One impor¬ 
tant point is that, e.g. unlike KMC, here detailed balance 
is not directly used to determine the transition prob¬ 
abilities, and in fact is in general not fulfilled. How¬ 
ever, detailed balance is naturally recovered at equilib¬ 
rium in calculations. Figure [5] shows example results for 
the equilibrium properties of the 2D Ising model with 
size = 100 X 100 using both ED and Metropolis, equi¬ 
librated for 5 X lO^F”^ and 5 x lO^MCS, respectively. 
The system undergoes a phase transition from param¬ 
agnetic to ferromagnetic phase at Tc = 2.269. The fig¬ 
ure shows that the algorithm reproduces well the results 
from Metropolis for the magnetization and the energy as 
a function of T, finding the same equilibrium configura¬ 
tions and Tc. 

While central to the algorithm, ted can also capture 
certain interesting aspects of the system dynamics. In 
the inset of Fig.[2l the characteristic time ted, averaged 
over a time at least equal to the equilibration one, is 
plotted as a function of T. For T > Tc, (ted) is very 
small ((ted) 10“^/F), corresponding to a fast flip- 
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ping rate, as expected in the paramagnetic phase. Below 
Tc, however, (ted) rapidly increases until it saturates 
for T < 1. At this temperature, the dynamics is essen¬ 
tially frozen and ted becomes one by construction. As 
discussed above, frozen dynamics is always reached in 
calculations for finite systems evolving into a stable con¬ 
figuration (e.g., the ferromagnetic state). This can often 
result in an unwanted slowing down of computations at 
sufficiently low T. A rapid growth of ted (e.g., below 
Tc in the figure) is then a computationally helpful flag 
of reaching a stable spin configuration. In fact, this is a 
limiting case of the time coarsening that is performed by 
Ted at each step of the algorithm, since ted is chosen to 
prevent spurious updating for At < ted at each step. 

We test the dynamical behavior of the algorithm for 
the 2D Ising model by quenching T from a disordered 
configuration (i.e., T = oo) to a subcritical temperature 
T = 1 < Tc corresponding to the fully magnetized ground 
state m Ki 1 (see Fig. [2]). This is a well known coarsen¬ 
ing process, where, as a result of quenching to low T, 
a mosaic of competing ordered-phase clusters is formed. 
In a finite system, the final state corresponds either to 
the fully ordered ground state or to a configuration with 
striped domains oriented antiparallel to the rest of the 
system [l^. The two physically relevant times in this 
situation are the time ti associated with the appearence 
of the first percolating cluster, i.e. an ordered domain 
of the size of the system, and the equilibration time Tgq 
after which the system is found in one of the two final 
states. 

The evolution of the magnetization m after the quench 
is shown in Fig. [3^, where results are averaged over 500 
quench realizations in systems of up to = 250 x 250 
spins. The hgure shows that the ED algorithm repro¬ 
duces the scaling m ~ typical of the coarsening dy¬ 
namics of two-dimensional systems with non-conserved 
order parameters Q, which is also captured by the 
Metropolis dynamics, reaching the equilibrium configu¬ 
ration (i.e., plateau in the figure) at Teq ~ 5 x 10^/F. 
The time t;, shown in the inset as a function of the lin¬ 
ear size T, signals the formation of the first percolating 
cluster, which has been demonstrated to be in general 
unstable 0- Its computation was performed by first 
determining the ferromagnetic clusters, using an imple¬ 
mentation of the Hoshen-Kopelman algorithm [l^, and 
then by checking the percolating properties along the ED 
dynamics. We find an exponent 6 = 0.84 for the power 
law Ti ~ T®, enriching the discussion on the phase order 
kinetics of models with non-conserved order parameter, 
usually developed within the KMC scheme [17|. 

Further information on the quench dynamics is ob¬ 
tained by the time evolution of (ted) shown in Fig. [SJd. 
Firstly, the equilibration time Teq extracted from Fig. [3^ 
is well captured by the dynamics of (ted). Consistently, 
the value (ted) ~ O.l/F for the plateau in Fig. [SJd is the 
same as that obtained at equilibrium for the correspond- 




FIG. 3: (a) Evolution of the magnetization m using ED (tri¬ 
angles) and Metropolis (squares), after a quench from a disor¬ 
dered configuration into T — 1, for a system of = 250x250. 
Inset: size scaling of the characteristic time ti at which the 
first percolation cluster is formed (see the text), (b) Evolu¬ 
tion of the average time step of the ED algorithm after the 
quench described in panel (a). Inset: size scaling of charac¬ 
teristic times r* and Teq represented in the main figure (see 
the text). 


ing temperature T = I (see inset of Fig. [2]). In addition, 
(ii) new information is provided by (ted) on the physical 
mechanisms of phase ordering. That is, a second char¬ 
acteristic time-scale t* appears at t* ^ 2 x lO^/F, just 
where (ted) changes the slope. By inspection, we find 
that T* corresponds to the appearance of the first few 
stationary (i.e., final) states in some realizations of the 
quenches. That is, no final configuration is reached in our 
simulations for t < t*. After this time, however, the sys¬ 
tem starts having a non-zero probability of being in the 
final state, where < ted > is maximal. Consequently, 
the average time scale of the relaxation slows down, in 
turn causing a more pronounced slope. In contrast, for 
t > Teq all configurations are either fully magnetized or 
striped. The inset shows that Teq (as well as t*) scales 
with the system size as Teq ^ T^, which is in agreement 
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with known results 112|. 


CONCLUSIONS 


In summary, we have introduced and tested a novel 
algorithm to simulate the stochastic dynamics of dis¬ 
crete variable models. To the best of our knowledge, 
this is the first Monte Carlo method involving the exact 
physical time at all scales, with no heuristic or phase- 
space assumptions. The latter opens up the study of, 
e.g., strongly out-of-equilibrium situations for which ex¬ 
act numerical calculations are currently not possible in 
short-time regimes. 

As said above, the present algorithm can be adapted to 
tackle several classes of different problems. For example, 
a microscopic update can be generalized that is consistent 
with conserved order parameter dynamics. The latter 
can describe, e.g., the dwiamics of kinetic phase separa¬ 
tion in binary mixtures [l[ . The role of two level subsys¬ 
tems is here played by each couple of nearest-neighbor 
sites with different occupations, while minimal processes 
translate into exchanges within these subsystems. The 
same reasoning applies to general Potts models and re¬ 
lated markovian chains. The study of quenches in classi¬ 
cal many-body systems and the relation to Kibble-Zurek 
mechanism |19l| . Lieb-Robinson bounds with short- and 
long-range interactions [ 2 ^, as well as dynamical phase 
transitions in magnetic models 2l|, are other important 
examples of physical processes of current interest where 
our algorithm can be straightforwardly applied. 
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